//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include "ephi.hpp"

CoilData::CoilData ()
{
  bmap = 0;
  bmapemap = 0;
  cmap = 0;
  cubic_bmap = 0;
  cubic_bmapemap = 0;
}

CoilData::~CoilData ()
{
  cleanup();
}

void CoilData::cleanup ()
{
  if ( bmap ) free (bmap);
  if ( bmapemap ) free (bmapemap);
  if ( cmap ) free (cmap);
  bmap = 0;
  bmapemap = 0;
  cmap = 0;
}

bool CoilData::load (const char *path)
{
  cleanup();

  FILE *fp = fopen (path, "rb");
  if ( !fp ) return false;

  char hdr_data[COILGEN_HEADER_SIZE];
  if ( !fread (hdr_data, sizeof (hdr_data), 1, fp) )
    {
      fclose (fp);
      return false;
    }

  const coilgen_header *hdr = (const coilgen_header*)hdr_data;

  if ( memcmp (hdr->magic, COILGEN_MAGIC, sizeof(hdr->magic)) != 0 )
      return false;

  flags = hdr->flags;

  coil_radius = hdr->coil_radius;
  coil_wr = hdr->coil_wr;
  coilpair_distance = hdr->coilpair_distance;

  if ( coil_radius < coil_wr || coil_radius < 0 || coil_wr < 0 )
      return false;
  
  radial_samples = hdr->radial_samples;
  axial_samples = hdr->axial_samples;
  geom_radial = hdr->geom_radial;
  geom_axial = hdr->geom_axial;
  geom_axial_add = hdr->geom_axial_add;

  has_axial_add = hdr->geom_axial_add != 0.0;

  mirror = (flags & coilgen_header::IS_MIRROR) ? true : false;

  cr = (prec_t)geom_radial / (prec_t)(radial_samples - COILGEN_BUFFER * 2 - 1);
  rec_cr = (prec_t)(radial_samples - COILGEN_BUFFER * 2 - 1) / (prec_t)geom_radial;
  ca = (prec_t)geom_axial / (prec_t)(axial_samples - COILGEN_BUFFER * 2 - 1);
  rec_ca = (prec_t)(axial_samples - COILGEN_BUFFER * 2 - 1) / (prec_t)geom_axial;

  Ephi::debug (Ephi::DEBUG_INFO, "Loading coil from %s, flags = %x\n", path, flags);

  while (1)
    {
      coilgen_blob_header bhdr;
      if ( !fread (&bhdr, sizeof(bhdr), 1, fp) )
        {
          // we assume this is the end of file
          fclose (fp);
          return true;
        }

      if ( !bhdr.can_read () )
        {
          Ephi::debug (Ephi::DEBUG_ERROR, "Unable to read blob, won't fit into memory\n");
          fclose (fp);
          return false;
        }

      switch (bhdr.blob_type)
        {
          case coilgen_blob_header::IS_BLOB_BMAP: {
            bmap = (coilgen_bmap *)malloc(bhdr.get_len());
            if ( !fread (bmap, bhdr.get_len(), 1, fp) )
              {
                Ephi::debug (Ephi::DEBUG_ERROR, "Failed to load bmap blob\n");
                fclose (fp);
                return false;
              }
            break;
          }
          case coilgen_blob_header::IS_BLOB_BMAPEMAP: {
            bmapemap = (coilgen_bmapemap *)malloc(bhdr.get_len());
            if ( !fread (bmapemap, bhdr.get_len(), 1, fp) )
              {
                Ephi::debug (Ephi::DEBUG_ERROR, "Failed to load emap blob\n");
                fclose (fp);
                return false;
              }
            break;
          }
          case coilgen_blob_header::IS_BLOB_CUBIC_BMAP: {
            cubic_bmap = (coilgen_cubic_bmap *)malloc(bhdr.get_len());
            if ( !fread (cubic_bmap, bhdr.get_len(), 1, fp) )
              {
                Ephi::debug (Ephi::DEBUG_ERROR, "Failed to load cubic bmap blob\n");
                fclose (fp);
                return false;
              }
            break;
          }
          case coilgen_blob_header::IS_BLOB_CUBIC_BMAPEMAP: {
            cubic_bmapemap = (coilgen_cubic_bmapemap *)malloc(bhdr.get_len());
            if ( !fread (cubic_bmapemap, bhdr.get_len(), 1, fp) )
              {
                Ephi::debug (Ephi::DEBUG_ERROR, "Failed to load cubic emap blob\n");
                fclose (fp);
                return false;
              }
            break;
          }
          case coilgen_blob_header::IS_BLOB_CONTACTMAP: {
            cmap = (unsigned char *)malloc(bhdr.get_len());
            if ( !fread (cmap, bhdr.get_len(), 1, fp) )
              {
                Ephi::debug (Ephi::DEBUG_ERROR, "Failed to load cmap blob of size %u\n", bhdr.get_len());
                fclose (fp);
                return false;
              }
            break;
          }
          default: {
            fseek (fp, bhdr.get_len(), SEEK_CUR);
            Ephi::debug (Ephi::DEBUG_WARN, "Unknown type blob %u found\n", bhdr.blob_type);
            break;
          }
        }
    }

  fclose (fp);

  return true;
}

// this is ugly

#define COIL_MATRIX_N 32
#define COIL_RADIUS 1.0
#define COIL_WR 0.1
#define MEMORY_USAGE 10 // megabytes
#define GEOM_RADIAL 2.0
#define GEOM_AXIAL 4.0

static size_t samplecount_radial;
static size_t samplecount_axial;
static prec_t cr, ca;
static prec_t geom_radial = GEOM_RADIAL;
static prec_t geom_axial = GEOM_AXIAL;
static prec_t geom_axial_add = 0.0;

static coilgen_bmap *bmap;
static coilgen_bmapemap *bmapemap;
static unsigned char *cmap;
static coilgen_cubic_bmap *cubic_bmap;
static coilgen_cubic_bmapemap *cubic_bmapemap;

struct coilgen_cubic_derivs
{
  struct {
    prec_t f;
    prec_t dfr;
    prec_t dfa;
    prec_t d2fra;
  } r, a;
};
static coilgen_cubic_derivs *bcubic;
static coilgen_cubic_derivs *ecubic;

static prec_t cubic_table[16][16] = {
  { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
  {-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0},
  { 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0},
  { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
  { 0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1},
  { 0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1},
  {-3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  { 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0},
  { 9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2},
  {-6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2},
  { 2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
  { 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0},
  {-6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1},
  { 4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1}
};

struct mytask : Task
{
  mytask(Statics& statics, ShapedCoilFactory::coils_t& coils, size_t si, size_t ei) :
      statics(statics), coils(coils), si(si), ei(ei) { }
  Statics& statics;
  ShapedCoilFactory::coils_t& coils;
  size_t si;
  size_t ei;

  void execute()
  {
    prec_t crc = cr * COILGEN_BUFFER;
    prec_t cac = ca * COILGEN_BUFFER;
    prec_t hr, ha;
    hr = ha = 5e-05;
    
    coilgen_bmap *bfw = 0;
    if ( bmap )
        bfw = bmap + si * samplecount_axial;
    coilgen_bmapemap *efw = 0;
    if ( bmapemap )
        efw = bmapemap + si * samplecount_axial;
    coilgen_cubic_derivs *bcd = 0;
    if ( bcubic )
        bcd = bcubic + si * samplecount_axial;
    coilgen_cubic_derivs *ecd = 0;
    if ( ecubic )
        ecd = ecubic + si * samplecount_axial;
    for ( size_t i = si; i < ei; i++ )
      {
        for ( size_t j = 0; j < samplecount_axial; j++ )
          {
            vect3d pos(i * cr - crc, 0, j * ca - cac - geom_axial_add);
            vect3d bf, ef;

            if ( bcd || ecd )
              {
//#define USE_NUMERICDIFF
#ifdef USE_NUMERICDIFF
                statics.getFields (pos, bf, ef);
                // we still interpolate the field value (consisting of two components:
                // radial and axial), so we have to calculate how both of these values
                // respond to changes in radial and axial _position_.
                prec_t rhr, rha;
                vect3d bf_rn, bf_rp, bf_an, bf_ap;
                vect3d ef_rn, ef_rp, ef_an, ef_ap;
                vect3d pa, pb;
                pa = vect3d(pos.x - hr, pos.y, pos.z);
                pb = vect3d(pos.x + hr, pos.y, pos.z);
                rhr = prec_t_sub (pb.x, pa.x);
                statics.getFields (pa, bf_rn, ef_rn);
                statics.getFields (pb, bf_rp, ef_rp);

                pa = vect3d(pos.x, pos.y, pos.z - ha);
                pb = vect3d(pos.x, pos.y, pos.z + ha);
                rha = prec_t_sub (pb.z, pa.z);
                statics.getFields (pa, bf_an, ef_an);
                statics.getFields (pb, bf_ap, ef_ap);

                // values for the cross-derivative, we call them in a ccw circle, 0, 1, 2, 3
                // where 0 is bottom-left corner
                vect3d bf0, bf1, bf2, bf3;
                vect3d ef0, ef1, ef2, ef3;
                pa = vect3d (pos.x - hr, pos.y, pos.z - ha);
                pb = vect3d (pos.x + hr, pos.y, pos.z - ha);
                statics.getFields (pa, bf0, ef0);
                statics.getFields (pb, bf1, ef1);
                pa.z = pb.z = pos.z + ha;
                statics.getFields (pa, bf3, ef3);
                statics.getFields (pb, bf2, ef2);
#endif

                if ( bcd )
                  {
#ifdef USE_NUMERICDIFF
                    bcd->r.f = prec2coilprec(bf.x);
                    bcd->r.dfr = (bf_rp.x - bf_rn.x) / rhr;
                    bcd->r.dfa = (bf_ap.x - bf_an.x) / rha;
                    bcd->r.d2fra = (bf2.x - bf3.x - bf1.x + bf0.x) / (rhr * rha);

                    bcd->a.f = prec2coilprec(bf.z);
                    bcd->a.dfr = (bf_rp.z - bf_rn.z) / rhr;
                    bcd->a.dfa = (bf_ap.z - bf_an.z) / rha;
                    bcd->a.d2fra = (bf2.z - bf3.z - bf1.z + bf0.z) / (rhr * rha);
#endif

#ifndef USE_NUMERICDIFF
                    prec_t cumBa, cumdBada, cumdBadr, cumd2Badadr, cumBr, cumdBrda, cumdBrdr, cumd2Brdadr;
                    cumBa = cumdBada = cumdBadr = cumd2Badadr = cumBr = cumdBrda = cumdBrdr = cumd2Brdadr = 0;
                    for ( size_t k = 0; k < coils.size(); k++ )
                      {
                        prec_t Ba, dBada, dBadr, d2Badadr, Br, dBrda, dBrdr, d2Brdadr;
                        calcCoilBField (coils[k].R, pos.z - coils[k].x, pos.x,
                            Ba, dBada, dBadr, d2Badadr, Br, dBrda, dBrdr, d2Brdadr);
                        cumBa += coils[k].I * Ba;
                        cumdBada += coils[k].I * dBada;
                        cumdBadr += coils[k].I * dBadr;
                        cumd2Badadr += coils[k].I * d2Badadr;
                        cumBr += coils[k].I * Br;
                        cumdBrda += coils[k].I * dBrda;
                        cumdBrdr += coils[k].I * dBrdr;
                        cumd2Brdadr += coils[k].I * d2Brdadr;
                      }

#if 0
                    char b[128];
                    printf ("%s: bf=(%g, %g), calc=(%g, %g)\n", pos.sprint(b), bf.z, bf.x, cumBa, cumBr);

                    printf ("%s: dada=%g, calc=%g\n", pos.sprint(b), bcd->a.dfa, cumdBada);
                    printf ("%s: dadr=%g, calc=%g\n", pos.sprint(b), bcd->a.dfr, cumdBadr);
                    printf ("%s: d2adadr=%g, calc=%g\n", pos.sprint(b), bcd->a.d2fra, cumd2Badadr);

                    printf ("%s: drda=%g, calc=%g\n", pos.sprint(b), bcd->r.dfa, cumdBrda);
                    printf ("%s: drdr=%g, calc=%g\n", pos.sprint(b), bcd->r.dfr, cumdBrdr);
                    printf ("%s: d2rdadr=%g, calc=%g\n", pos.sprint(b), bcd->r.d2fra, cumd2Brdadr);
#endif

                    bcd->a.f = cumBa;
                    bcd->a.dfa = cumdBada;
                    bcd->a.dfr = cumdBadr;
                    bcd->a.d2fra = cumd2Badadr;

                    bcd->r.f = cumBr;
                    bcd->r.dfa = cumdBrda;
                    bcd->r.dfr = cumdBrdr;
                    bcd->r.d2fra = cumd2Brdadr;
#endif

                    bcd++;
                  }

                if ( ecd )
                  {
#ifdef USE_NUMERICDIFF
                    ecd->r.f = prec2coilprec(ef.x);
                    ecd->r.dfr = (ef_rp.x - ef_rn.x) / rhr;
                    ecd->r.dfa = (ef_ap.x - ef_an.x) / rha;
                    ecd->r.d2fra = (ef2.x - ef3.x - ef1.x + ef0.x) / (rhr * rha);

                    ecd->a.f = prec2coilprec(ef.z);
                    ecd->a.dfr = (ef_rp.z - ef_rn.z) / rhr;
                    ecd->a.dfa = (ef_ap.z - ef_an.z) / rha;
                    ecd->a.d2fra = (ef2.z - ef3.z - ef1.z + ef0.z) / (rhr * rha);
#endif
#ifndef USE_NUMERICDIFF
                    prec_t cumEa, cumdEada, cumdEadr, cumd2Eadadr, cumEr, cumdErda, cumdErdr, cumd2Erdadr;
                    cumEa = cumdEada = cumdEadr = cumd2Eadadr = cumEr = cumdErda = cumdErdr = cumd2Erdadr = 0;
                    for ( size_t k = 0; k < coils.size(); k++ )
                      {
                        prec_t Ea, dEada, dEadr, d2Eadadr, Er, dErda, dErdr, d2Erdadr;
                        calcCoilEField (coils[k].R, pos.z - coils[k].x, pos.x,
                            Ea, dEada, dEadr, d2Eadadr, Er, dErda, dErdr, d2Erdadr);
                        cumEa += coils[k].q * Ea;
                        cumdEada += coils[k].q * dEada;
                        cumdEadr += coils[k].q * dEadr;
                        cumd2Eadadr += coils[k].q * d2Eadadr;
                        cumEr += coils[k].q * Er;
                        cumdErda += coils[k].q * dErda;
                        cumdErdr += coils[k].q * dErdr;
                        cumd2Erdadr += coils[k].q * d2Erdadr;
                      }

#if 0
                    char b[128];
                    printf ("%s: bf=(%g, %g), calc=(%g, %g)\n", pos.sprint(b), ef.z, ef.x, cumEa, cumEr);

                    printf ("%s: dada=%g, calc=%g\n", pos.sprint(b), bcd->a.dfa, cumdEada);
                    printf ("%s: dadr=%g, calc=%g\n", pos.sprint(b), bcd->a.dfr, cumdEadr);
                    printf ("%s: d2adadr=%g, calc=%g\n", pos.sprint(b), bcd->a.d2fra, cumd2Eadadr);

                    printf ("%s: drda=%g, calc=%g\n", pos.sprint(b), bcd->r.dfa, cumdErda);
                    printf ("%s: drdr=%g, calc=%g\n", pos.sprint(b), bcd->r.dfr, cumdErdr);
                    printf ("%s: d2rdadr=%g, calc=%g\n", pos.sprint(b), bcd->r.d2fra, cumd2Erdadr);
#endif

                    ecd->a.f = cumEa;
                    ecd->a.dfa = cumdEada;
                    ecd->a.dfr = cumdEadr;
                    ecd->a.d2fra = cumd2Eadadr;

                    ecd->r.f = cumEr;
                    ecd->r.dfa = cumdErda;
                    ecd->r.dfr = cumdErdr;
                    ecd->r.d2fra = cumd2Erdadr;
#endif
                    ecd++;
                  }
              }
            else
              {
                statics.getFields (pos, bf, ef);
                if ( bfw )
                  {
                    bfw->rb = prec2coilprec(bf.x);
                    bfw->ab = prec2coilprec(bf.z);
                    bfw++;
                  }
                else if ( efw )
                  {
                    efw->rb = prec2coilprec(bf.x);
                    efw->ab = prec2coilprec(bf.z);
                    efw->re = prec2coilprec(ef.x);
                    efw->ae = prec2coilprec(ef.z);
                    efw++;
                  }
              }
//            pos.print();
//            printf ("%lu: %g, %g\n", (unsigned long)(i * samplecount_axial + j), bfw->radial, bfw->axial);
            if ( statics.inContact (pos) )
              {
                size_t index = i * samplecount_axial + j;
//                printf ("%lu: %g, %g, in contact\n", (unsigned long)index, bfw->radial, bfw->axial);
                cmap[index >> 3] |= 1 << (index & 7);
              }
          }
      }
  }
};

struct mytask3 : Task
{
  mytask3(size_t si, size_t ei) : si(si), ei(ei) { }
  size_t si;
  size_t ei;

  void execute()
  {
    coilgen_cubic_derivs *bcd = 0;
    if ( bcubic )
        bcd = bcubic + si * samplecount_axial;
    coilgen_cubic_derivs *ecd = 0;
    if ( ecubic )
        ecd = ecubic + si * samplecount_axial;

    coilgen_cubic_bmap *bfw = 0;
    if ( cubic_bmap )
        bfw = cubic_bmap + si * (samplecount_axial - 1);
    coilgen_cubic_bmapemap *efw = 0;
    if ( cubic_bmapemap )
        efw = cubic_bmapemap + si * (samplecount_axial - 1);

    prec_t cacr = ca * cr;

    for ( size_t i = si; i < ei; i++ )
      {
        for ( size_t j = 0; j < samplecount_axial - 1; j++ )
          {
            coilgen_cubic_derivs *b0 = bcd;
            coilgen_cubic_derivs *b1 = bcd + 1;
            coilgen_cubic_derivs *b3 = bcd + samplecount_axial;
            coilgen_cubic_derivs *b2 = bcd + samplecount_axial + 1;
            coilgen_cubic_derivs *ba[4] = {b0, b1, b2, b3};

#if 0
            printf ("b0(%p): %g, %g, %g, %g; %g, %g, %g, %g\n",
                b0,
                b0->r.f, b0->r.dfr, b0->r.dfa, b0->r.d2fra,
                b0->a.f, b0->a.dfr, b0->a.dfa, b0->a.d2fra);
#endif

            prec_t v[16];
            prec_t r[16];
            for ( size_t i = 0; i < 4; i++ )
              {
                v[i] = ba[i]->r.f;
                v[i+4] = ba[i]->r.dfa * ca;
                v[i+8] = ba[i]->r.dfr * cr;
                v[i+12] = ba[i]->r.d2fra * cacr;
              }
            for ( size_t i = 0; i < 16; i++ )
              {
                prec_t cum = 0;
                for ( int j = 0; j < 16; j++ )
                    cum += v[j] * cubic_table[i][j];
                r[i] = cum;
              }
            for ( size_t i = 0; i < 4; i++ )
              {
                for ( size_t j = 0; j < 4; j++ )
                  {
                    if ( bfw )
                        bfw->rb[i][j] = r[i * 4 + j];
                    if ( efw )
                        efw->rb[i][j] = r[i * 4 + j];
                  }
              }

            for ( size_t i = 0; i < 4; i++ )
              {
                v[i] = ba[i]->a.f;
                v[i+4] = ba[i]->a.dfa * ca;
                v[i+8] = ba[i]->a.dfr * cr;
                v[i+12] = ba[i]->a.d2fra * cacr;
              }
            for ( size_t i = 0; i < 16; i++ )
              {
                prec_t cum = 0;
                for ( int j = 0; j < 16; j++ )
                    cum += v[j] * cubic_table[i][j];
                r[i] = cum;
              }
            for ( size_t i = 0; i < 4; i++ )
              {
                for ( size_t j = 0; j < 4; j++ )
                  {
                    if ( bfw )
                        bfw->ab[i][j] = r[i * 4 + j];
                    if ( efw )
                        efw->ab[i][j] = r[i * 4 + j];
                  }
              }

            if ( efw )
              {
                coilgen_cubic_derivs *e0 = ecd;
                coilgen_cubic_derivs *e1 = ecd + 1;
                coilgen_cubic_derivs *e3 = ecd + samplecount_axial;
                coilgen_cubic_derivs *e2 = ecd + samplecount_axial + 1;
                coilgen_cubic_derivs *ea[4] = {e0, e1, e2, e3};

                for ( size_t i = 0; i < 4; i++ )
                  {
                    v[i] = ea[i]->r.f;
                    v[i+4] = ea[i]->r.dfa * ca;
                    v[i+8] = ea[i]->r.dfr * cr;
                    v[i+12] = ea[i]->r.d2fra * cacr;
                  }
                for ( size_t i = 0; i < 16; i++ )
                  {
                    prec_t cum = 0;
                    for ( int j = 0; j < 16; j++ )
                        cum += v[j] * cubic_table[i][j];
                    r[i] = cum;
                  }
                for ( size_t i = 0; i < 4; i++ )
                    for ( size_t j = 0; j < 4; j++ )
                        efw->re[i][j] = r[i * 4 + j];

                for ( size_t i = 0; i < 4; i++ )
                  {
                    v[i] = ea[i]->a.f;
                    v[i+4] = ea[i]->a.dfa * ca;
                    v[i+8] = ea[i]->a.dfr * cr;
                    v[i+12] = ea[i]->a.d2fra * cacr;
                  }
                for ( size_t i = 0; i < 16; i++ )
                  {
                    prec_t cum = 0;
                    for ( int j = 0; j < 16; j++ )
                        cum += v[j] * cubic_table[i][j];
                    r[i] = cum;
                  }
                for ( size_t i = 0; i < 4; i++ )
                    for ( size_t j = 0; j < 4; j++ )
                        efw->ae[i][j] = r[i * 4 + j];
                efw++;
              }
            if ( bfw ) bfw++;
            bcd++;
            ecd++;
          }
      }
  }
};

CoilGen::CoilGen ()
{
  geom_radial = GEOM_RADIAL;
  geom_axial = GEOM_AXIAL;
  radius = COIL_RADIUS;

  wr = COIL_WR;
  dbl = 0;
  shape = ROUND;
  memory_usage = MEMORY_USAGE;
  coil_matrix_n = COIL_MATRIX_N;
  do_efield = false;
  do_cubic = false;

  force_grid = 0;
}

void CoilGen::generate(const char *path, bool print_info, bool show_progress)
{
  Statics statics;

  ::geom_radial = geom_radial;
  ::geom_axial = geom_axial;
  geom_axial_add = 0;

  size_t entrysize;
  if ( do_cubic )
    {
      entrysize = do_efield ? sizeof (coilgen_cubic_bmapemap) : sizeof (coilgen_cubic_bmap);
    }
  else
    {
      entrysize = do_efield ? sizeof (coilgen_bmapemap) : sizeof (coilgen_bmap);
    }
  size_t mem = memory_usage * 1024 * 1024;
  size_t samples = mem/entrysize;
  samplecount_radial = prec2size_t(prec_t_sqrt ((prec_t)samples * (geom_radial / geom_axial)));
  samplecount_axial = prec2size_t(prec_t_sqrt ((prec_t)samples / (geom_radial / geom_axial)));
  if ( force_grid )
    {
      samplecount_radial = force_grid;
      samplecount_axial = force_grid;
    }
  cr = (prec_t)geom_radial / (prec_t)(samplecount_radial - 1);
  ca = (prec_t)geom_axial / (prec_t)(samplecount_axial - 1);
  samplecount_radial += 2 * COILGEN_BUFFER;
  samplecount_axial += 2 * COILGEN_BUFFER;

  size_t fmap_size;
  if ( do_cubic )
    {
      fmap_size = entrysize * (samplecount_radial - 1) * (samplecount_axial - 1);
      if ( do_efield )
        {
          ecubic = new coilgen_cubic_derivs [samplecount_radial * samplecount_axial];
          cubic_bmapemap = new coilgen_cubic_bmapemap [(samplecount_radial - 1) * (samplecount_axial - 1)];
        }
      else
        {
          cubic_bmap = new coilgen_cubic_bmap [(samplecount_radial - 1) * (samplecount_axial - 1)];
        }
      bcubic = new coilgen_cubic_derivs [samplecount_radial * samplecount_axial];
    }
  else
    {
      fmap_size = entrysize * samplecount_radial * samplecount_axial;
      if ( do_efield )
        {
          bmap = 0;
          bmapemap = new coilgen_bmapemap [samplecount_radial * samplecount_axial];
        }
      else
        {
          bmap = new coilgen_bmap [samplecount_radial * samplecount_axial];
          bmapemap = 0;
        }
    }
  int cmap_size = samplecount_radial * samplecount_axial / 8 + 1;
  cmap = new unsigned char [cmap_size];
  memset (cmap, 0, samplecount_radial * samplecount_axial / 8 + 1);

  ShapedCoilFactory cf(coil_matrix_n);
  if ( coil_matrix_n > 1 )
    {
      if ( shape == ROUND )
          cf.set_circle();
      else
          cf.set_scale (1.0 / PREC_SQRT2, 1.0 / PREC_SQRT2);
    }

  ShapedCoilFactory::coils_t coils;
  prec_t q = do_efield ? COILGEN_QMOD : 0;
  if ( dbl != 0 )
    {
      vect3d target(0, 0, 0);
      if ( dbl < 0 )
          target = vect3d (0, 0, dbl);
      vect3d pos;
      
      pos = vect3d (0, 0, fabs(dbl) * 0.5);
      cf.make (statics, pos, target - pos, radius, wr, 1.0, q);
      cf.add (coils, (target - pos).z, radius, wr, 1.0, q);
      pos = vect3d (0, 0, -fabs(dbl) * 0.5);
      cf.make (statics, pos, target - pos, radius, wr, 1.0, q);
      cf.add (coils, (target - pos).z, radius, wr, -1.0, q);

      //geom_axial_add = geom_axial * 0.5;
    }
  else
    {
      cf.make (statics, vect3d(0, 0, 0), vect3d(0, 0, 1), radius, wr, 1.0, q);
      cf.add (coils, 0, radius, wr, 1.0, q);
    }

  if ( print_info )
    {
      printf ("Coil radius: %f m\n", prec2double (radius));
      printf ("Coil cross-section radius: %f m\n", prec2double(wr));
      switch (shape)
        {
          case ROUND: printf ("Coil shape: ROUND\n"); break;
          case SQUARE: printf ("Coil shape: SQUARE\n"); break;
        }
      printf ("Coil memoy allocation: %lu MB\n", (unsigned long)memory_usage);
      printf ("Plane square calculation geometry (radial x axial): %f x %f m\n",
          prec2double(geom_radial), prec2double(geom_axial));
      printf ("Coil matrix formed out of %lu x %lu coils\n", (unsigned long)coil_matrix_n, (unsigned long)coil_matrix_n);
      if ( dbl != 0 )
          printf ("Double coils, oriented %s, apart by %f m\n", dbl < 0 ? "same way" : "opposing", prec2double(fabs(dbl)));
      printf ("Sample field: %lu x %lu samples\n", (unsigned long)samplecount_radial, (unsigned long)samplecount_axial);
      printf ("Sample cell size: %g x %g m\n", prec2double(cr), prec2double(ca));
      printf ("Calculating the e-field: %s\n", do_efield ? "YES" : "NO");
      printf ("Calcule for cubic interpolation: %s\n", do_cubic ? "YES" : "NO");
      printf ("Writing data into %s\n", path);
      printf ("geom_axial_add = %f\n", prec2double(geom_axial_add));
    }

  TaskManager tmgr;
  for ( size_t i = 0; i < samplecount_radial; i++ )
      tmgr.addTask (new mytask(statics, coils, i, i + 1));
  size_t n = 0;
  if ( show_progress )
    {
      printf ("Calculating the grid..");
      fflush (stdout);
    }
  while (tmgr.run())
    {
      n++;
      if ( show_progress )
        {
          printf ("\rCalculating the grid: %lu%%", (unsigned long)(100 * n / samplecount_radial));
          fflush (stdout);
        }
    }
  if ( show_progress )
      printf ("\rGrid done!                           \n");

  if ( do_cubic )
    {
      for ( size_t i = 0; i < samplecount_radial - 1; i++ )
          tmgr.addTask (new mytask3(i, i + 1));
      size_t n = 0;
      if ( show_progress )
        {
          printf ("Calculating the cubic interpolation coefficents..");
          fflush (stdout);
        }
      while (tmgr.run())
        {
          n++;
          if ( show_progress )
            {
              printf ("\rCalculating the cubic interpolation coefficents: %lu%%",
                  (unsigned long)(100 * n / samplecount_radial));
              fflush (stdout);
            }
        }
      if ( show_progress )
          printf ("\rCubic done!                                                    \n");
    }

  char header[COILGEN_HEADER_SIZE];
  memset (header, 0, sizeof(header));
  coilgen_header *hdr = (coilgen_header*)header;
  memcpy (hdr->magic, COILGEN_MAGIC, sizeof (hdr->magic));
  hdr->flags = 0;
  if ( dbl < 0 )
      hdr->flags |= coilgen_header::IS_COILPAIR;
  else if ( dbl > 0 )
      hdr->flags |= coilgen_header::IS_COILPAIR | coilgen_header::IS_OPPOSING_COILPAIR;
  if ( geom_axial_add == 0.0 )
    {
      hdr->flags |= coilgen_header::IS_HALF;
      if ( dbl > 0 )
          hdr->flags |= coilgen_header::IS_MIRROR;
    }

  hdr->coilprec_size = sizeof (coilprec_t);

  hdr->coil_radius = prec2double(radius);
  hdr->coil_wr = prec2double(wr);
  hdr->coilpair_distance = prec2double(fabs(dbl));

  hdr->radial_samples = samplecount_radial;
  hdr->axial_samples = samplecount_axial;
  hdr->geom_radial = prec2double(geom_radial);
  hdr->geom_axial = prec2double(geom_axial);
  hdr->geom_axial_add = prec2double(geom_axial_add);

  FILE *fp = fopen (path, "wb");
  if ( !fp )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "Unable to open %s for writing!\n", path);
      exit(-1);
    }

  if ( !fwrite (header, sizeof (header), 1, fp) )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "Unable to write into %s!\n", path);
      exit(-1);
    }

  struct coilgen_blob_header bhdr;
  memset (&bhdr, 0, sizeof(bhdr));

  if ( do_efield )
    {
      if ( do_cubic )
        {
          bhdr.blob_type = coilgen_blob_header::IS_BLOB_CUBIC_BMAPEMAP;
          bhdr.set_len (fmap_size);
          fwrite (&bhdr, sizeof (bhdr), 1, fp);
          fwrite (cubic_bmapemap, fmap_size, 1, fp);
        }
      else
        {
          bhdr.blob_type = coilgen_blob_header::IS_BLOB_BMAPEMAP;
          bhdr.set_len (fmap_size);
          fwrite (&bhdr, sizeof (bhdr), 1, fp);
          fwrite (bmapemap, entrysize * samplecount_radial * samplecount_axial, 1, fp);
        }
    }
  else
    {
      if ( do_cubic )
        {
          bhdr.blob_type = coilgen_blob_header::IS_BLOB_CUBIC_BMAP;
          bhdr.set_len (fmap_size);
          fwrite (&bhdr, sizeof (bhdr), 1, fp);
          fwrite (cubic_bmap, fmap_size, 1, fp);
        }
      else
        {
          bhdr.blob_type = coilgen_blob_header::IS_BLOB_BMAP;
          bhdr.set_len (fmap_size);
          fwrite (&bhdr, sizeof (bhdr), 1, fp);
          fwrite (bmap, entrysize * samplecount_radial * samplecount_axial, 1, fp);
        }
    }

  bhdr.blob_type = coilgen_blob_header::IS_BLOB_CONTACTMAP;
  bhdr.set_len (cmap_size);
  fwrite (&bhdr, sizeof (bhdr), 1, fp);
  fwrite (cmap, cmap_size, 1, fp);
  fclose (fp);

  delete [] bmap;
  delete [] bmapemap;
  delete [] cubic_bmap;
  delete [] cubic_bmapemap;
  delete [] cmap;
  delete [] ecubic;
  delete [] bcubic;
  bmap = 0;
  bmapemap = 0;
  cubic_bmap = 0;
  cubic_bmapemap = 0;
  cmap = 0;
  ecubic = 0;
  bcubic = 0;
}

